1 Data

To run this vignette, we will use the Skin dataset from the article High-plex protein and whole transcriptome co-mapping at cellular resolution with spatial CITE-seq

2 Installation

To run this vignette, you must download the last version of Giotto Suite

remotes::install_github("drieslab/Giotto@suite_dev")

3 Individual modalities preparation

3.1 Create Giotto object

Load Giotto

library(Giotto)
## Giotto Suite 3.3.1
instrs <- createGiottoInstructions(save_plot = FALSE,
                                   show_plot = FALSE)

Create Giotto object using RNA and Protein expression, as well as spatial positions.

my_giotto_object <- createGiottoObject(
  expression = list(rna = list(raw = "data/expression_rna.csv"),
                    protein = list(raw = "data/expression_protein.csv")),
  expression_feat = list("rna", "protein"),
  spatial_locs = "data/positions.csv",
  instructions = instrs)

3.2 Add tissue image

my_giotto_image <- createGiottoImage(gobject = my_giotto_object,
                                     do_manual_adj = TRUE,
                                     scale_factor = 0.0625,
                                     mg_object = "data/skin.jpg",
                                     negative_y = TRUE)
## do_manual_adj == TRUE 
##  Boundaries will be adjusted by given values.
my_giotto_object <- addGiottoImage(gobject = my_giotto_object,
                                   images = list(my_giotto_image),
                                   spat_loc_name = "raw")

Visualize image

spatPlot2D(my_giotto_object,
         show_image = TRUE,
         point_size = 2)

3.3 RNA modality

3.3.1 Filtering

my_giotto_object <- filterGiotto(gobject = my_giotto_object,
                                 spat_unit = "cell",
                                 feat_type = "rna",
                                 expression_threshold = 1,
                                 feat_det_in_min_cells = 1,
                                 min_det_feats_per_cell = 1)
## completed 1: preparation 
## completed 2: subset expression data 
## completed 3: subset spatial locations 
## completed 4: subset cell (spatial units) and feature IDs 
## completed 5: subset cell metadata 
## completed 6: subset feature metadata 
## completed 7: subset spatial network(s) 
## completed 8: subsetted dimension reductions 
## completed 9: subsetted nearest network(s) 
## completed 10: subsetted spatial enrichment results 
## 
## Feature type:  rna 
## Number of cells removed:  0  out of  1691 
## Number of feats removed:  363  out of  15486

3.3.2 Normalization

my_giotto_object <- normalizeGiotto(gobject = my_giotto_object,
                                    spat_unit = "cell",
                                    feat_type = "rna",
                                    norm_methods = "standard",
                                    scalefactor = 10000,
                                    verbose = TRUE)
## 
## first scale feats and then cells

3.3.3 Gene statistics

my_giotto_object <- addStatistics(gobject = my_giotto_object,
                                  spat_unit = "cell",
                                  feat_type = "rna")

3.3.4 HVF identification

my_giotto_object <- calculateHVF(gobject = my_giotto_object,
                                 spat_unit = "cell",
                                 feat_type = "rna",
                                 expression_values = "normalized",
                                 method = "cov_groups",
                                 nr_expression_groups = 20,
                                 zscore_threshold = 1.5)
## return_plot = TRUE and return_gobject = TRUE 
## 
##           plot will not be returned to object, but can still be saved with save_plot = TRUE or manually

3.3.5 Principal component analysis (PCA)

my_giotto_object <- runPCA(gobject = my_giotto_object,
                           spat_unit = "cell",
                           feat_type = "rna",
                           expression_values = "normalized",
                           reduction = "cells",
                           feats_to_use = "hvf",
                           name = "rna.pca")
## "hvf" was found in the feats metadata information and will be used to select
##  highly variable features
## class of selected matrix:  dgCMatrix 
## [1] "finished runPCA_factominer, method == factominer"

3.3.6 Uniform manifold approximation projection (UMAP)

my_giotto_object <- runUMAP(gobject = my_giotto_object,
                            spat_unit = "cell",
                            feat_type = "rna",
                            expression_values = "normalized",
                            reduction = "cells",
                            dimensions_to_use = 1:10,
                            dim_reduction_name = "rna.pca")

3.3.7 Create nearest network

my_giotto_object <- createNearestNetwork(gobject = my_giotto_object,
                                         spat_unit = "cell",
                                         feat_type = "rna",
                                         type = "kNN",
                                         dim_reduction_name = "rna.pca",
                                         name = "rna_kNN.pca",
                                         dimensions_to_use = 1:10,
                                         k = 20)

3.3.8 Find Leiden clusters

my_giotto_object <- doLeidenCluster(gobject = my_giotto_object,
                                    spat_unit = "cell",
                                    feat_type = "rna",
                                    nn_network_to_use = "kNN",
                                    network_name = "rna_kNN.pca",
                                    name = "leiden_clus")

3.3.9 Plot PCA

plotPCA(gobject = my_giotto_object,
        spat_unit = "cell",
        feat_type = "rna",
        dim_reduction_name = "rna.pca",
        cell_color = 'leiden_clus',
        title = "RNA PCA")

3.3.10 Plot UMAP

plotUMAP(gobject = my_giotto_object,
         spat_unit = "cell",
         feat_type = "rna",
         cell_color = 'leiden_clus',
         point_size = 2,
         title = "RNA Uniform Manifold Approximation & Projection (UMAP)",
         axis_title = 12,
         axis_text = 10)

3.3.11 Plot spatial locations by cluster

spatPlot2D(my_giotto_object,
           show_image = TRUE,
           point_size = 1.5,
           cell_color = "leiden_clus",
           title = "RNA Leiden clustering")

3.4 Protein modality

3.4.1 Filtering

my_giotto_object <- filterGiotto(gobject = my_giotto_object,
                                 spat_unit = "cell",
                                 feat_type = "protein",
                                 expression_threshold = 1,
                                 feat_det_in_min_cells = 1,
                                 min_det_feats_per_cell = 1)
## completed 1: preparation 
## completed 2: subset expression data 
## completed 3: subset spatial locations 
## completed 4: subset cell (spatial units) and feature IDs 
## completed 5: subset cell metadata 
## completed 6: subset feature metadata 
## completed 7: subset spatial network(s) 
## completed 8: subsetted dimension reductions 
## completed 9: subsetted nearest network(s) 
## completed 10: subsetted spatial enrichment results 
## 
## Feature type:  protein 
## Number of cells removed:  0  out of  1691 
## Number of feats removed:  0  out of  283

3.4.2 Normalization

my_giotto_object <- normalizeGiotto(gobject = my_giotto_object,
                                    spat_unit = "cell",
                                    feat_type = "protein",
                                    scalefactor = 10000,
                                    verbose = T)
## 
## first scale feats and then cells

3.4.3 Gene statistics

my_giotto_object <- addStatistics(gobject = my_giotto_object,
                                  spat_unit = "cell",
                                  feat_type = "protein",
                                  expression_values = "normalized")

3.4.4 Principal component analysis (PCA)

my_giotto_object <- runPCA(gobject = my_giotto_object,
                           spat_unit = "cell",
                           feat_type = "protein",
                           expression_values = "normalized",
                           scale_unit = T,
                           center = F,
                           method = "factominer")
## "hvf" was not found in the gene metadata information, all genes will be used
## class of selected matrix:  dgCMatrix 
## [1] "finished runPCA_factominer, method == factominer"

3.4.5 Uniform manifold approximation projection (UMAP)

my_giotto_object <- runUMAP(gobject = my_giotto_object,
                            spat_unit = "cell",
                            feat_type = "protein",
                            expression_values = "normalized",
                            dimensions_to_use = 1:10)

3.4.6 Create nearest network

my_giotto_object <- createNearestNetwork(gobject = my_giotto_object,
                                         spat_unit = "cell",
                                         feat_type = "protein",
                                         type = "kNN",
                                         name = "protein_kNN.pca",
                                         dimensions_to_use = 1:10,
                                         k = 20)

3.4.7 Find Leiden clusters

my_giotto_object <- doLeidenCluster(gobject = my_giotto_object,
                                    spat_unit = "cell",
                                    feat_type = "protein",
                                    nn_network_to_use = "kNN",
                                    network_name = "protein_kNN.pca",
                                    name = "leiden_clus")

3.4.8 Plot PCA

plotPCA(gobject = my_giotto_object,
        spat_unit = "cell",
        feat_type = "protein",
        dim_reduction_name = "protein.pca",
        cell_color = 'leiden_clus',
        title = "Protein PCA")

3.4.9 Plot UMAP

plotUMAP(gobject = my_giotto_object,
         spat_unit = "cell",
         feat_type = "protein",
         cell_color = 'leiden_clus',
         dim_reduction_name = "protein.umap",
         point_size = 2,
         title = "Protein Uniform Manifold Approximation & Projection (UMAP)",
         axis_title = 12,
         axis_text = 10)

3.4.10 Plot spatial locations by cluster

spatPlot2D(my_giotto_object,
           spat_unit = "cell",
           feat_type = "protein",
           cell_color = "leiden_clus",
           point_size = 1.5,
           show_image = TRUE,
           title = "Protein Leiden clustering")

3.5 Save Giotto object

saveGiotto(my_giotto_object, "multiomics_Giotto_object")

4 Multi-omics integration

4.1 Load Giotto object

Download the file from https://drive.google.com/file/d/1GyhKWl14MxEJCw_2-e8Bge3d3zY7XLkc/view?usp=share_link

my_giotto_object <- loadGiotto("multiomics_Giotto_object")

4.2 Calculate WNN

my_giotto_object <- runWNN(my_giotto_object,
                           modality_1 = "rna",
                           modality_2 = "protein",
                           pca_name_modality_1 = "rna.pca",
                           pca_name_modality_2 = "protein.pca",
                           k = 20,
                           verbose = TRUE)
## The NN network name was not specified, default to the first: "rna_kNN.pca"
## The NN network name was not specified, default to the first: "protein_kNN.pca"
## [1] "Calculating rna - rna distance"
## [1] "Calculating protein - protein distance"
## [1] "Calculating low dimensional cell-cell distances for rna"
## [1] "Calculating low dimensional cell-cell distances for protein"
## [1] "Calculating within-modality prediction"
## [1] "Calculating cross-modality prediction"
## [1] "Calculating Jaccard similarities"
## [1] "Calculating kernel bandwidths"
## [1] "Calculating modality weights"
## [1] "Calculating WNN"
## [1] "Calculating WNN normalization"
## [1] "Calculating WNN graph"
## [1] "Saving WNN results"

4.3 Create integrated UMAP

my_giotto_object <- runIntegratedUMAP(my_giotto_object,
                                      modality1 = "rna",
                                      modality2 = "protein",
                                      spread = 5,
                                      min_dist = 2)
## [1] "Calculating integrated Nearest Neighbors"
## [1] "Creating igraph"
## [1] "Calculating integrated UMAP"

4.4 Calculate Leiden clusters

my_giotto_object <- doLeidenCluster(gobject = my_giotto_object,
                                    spat_unit = "cell",
                                    feat_type = "rna",
                                    nn_network_to_use = "kNN",
                                    network_name = "integrated_kNN",
                                    name = "integrated_leiden_clus",
                                    resolution = 0.4)

4.5 Plot integrated UMAP

plotUMAP(gobject = my_giotto_object,
         spat_unit = "cell",
         feat_type = "rna",
         cell_color = 'integrated_leiden_clus',
         dim_reduction_name = "integrated.umap",
         point_size = 1.5,
         title = "Integrated UMAP using Integrated Leiden clusters",
         axis_title = 12,
         axis_text = 10)

4.6 Plot integrated spatial locations by cluster

spatPlot2D(my_giotto_object,
           spat_unit = "cell",
           feat_type = "rna",
           cell_color = "integrated_leiden_clus",
           point_size = 1.5,
           show_image = TRUE,
           title = "Integrated Leiden clustering")

5 Interactive selection

5.1 Add tissue image

my_giotto_image <- createGiottoImage(gobject = my_giotto_object,
                                     do_manual_adj = TRUE,
                                     scale_factor = 0.0625,
                                     mg_object = "data/skin.jpg",
                                     negative_y = TRUE)

my_giotto_object <- addGiottoImage(gobject = my_giotto_object,
                                   images = list(my_giotto_image),
                                   spat_loc_name = "raw")

The interactive selection tool will allow you to run a Shiny app and plot polygons over the tissue.

Note that you can use the blue bars to Zoom in/out at areas of interest over the tissue

You can choose multiple polygon names, one for each drawn polygon.

When you have finished plotting polygons, click on the Done button

my_polygon_coordinates will store a data.table object with the polygon names and xy coordinates for downstream analysis.

5.2 Create a spatial plot

my_spatPlot <- spatPlot2D(gobject = my_giotto_object,
                          spat_unit = "cell",
                          feat_type = "protein",
                          show_image = TRUE,
                          cell_color = 'leiden_clus',
                          point_size = 2,
                          point_alpha = 0.75)

my_spatPlot

5.3 Run interactive tool

my_polygon_coordinates <- plotInteractivePolygons(my_spatPlot, color = "black")

5.4 Add polygons information to Giotto object

my_giotto_polygons <- createGiottoPolygonsFromDfr(my_polygon_coordinates, 
                                                  name = 'selections')
##   Selecting col "name" as poly_ID column
##   Selecting cols "x" and "y" as x and y respectively
my_giotto_object <- addGiottoPolygons(gobject = my_giotto_object, 
                                      gpolygons = list(my_giotto_polygons))

addPolygonCells will modify the cell_metadata slot to indicate if a certain cell is located within a polygon area (polygon1,2,3..) or not (no_polygon)

my_giotto_object <- addPolygonCells(my_giotto_object, 
                                    polygon_name = 'selections')
## 
## These column names were already used: nr_feats perc_feats total_expr
##  leiden_clus integrated_leiden_clus
##  and will be overwritten

Verify that cell metadata has been modified

head(pDataDT(my_giotto_object))
##    cell_ID nr_feats perc_feats total_expr leiden_clus integrated_leiden_clus
## 1:   10x14       70  0.4628711   474.8438           2                      7
## 2:   10x15      237  1.5671494  1228.7768           4                      7
## 3:   10x16       92  0.6083449   599.3736           2                      5
## 4:   10x17       87  0.5752827   581.0906           5                      4
## 5:   10x18      132  0.8728427   807.0911           7                      1
## 6:   10x19      336  2.2217814  1601.3652           3                      1
##    selections
## 1:  polygon 2
## 2:  polygon 2
## 3:  polygon 2
## 4:  polygon 2
## 5: no_polygon
## 6: no_polygon

5.5 Retrieve cells from the polygon area

my_polygons_cells is a terra::spatVector object that contains the information of each cell within the polygon areas

my_polygons_cells <- getCellsFromPolygon(my_giotto_object, 
                                         polygon_name = 'selections')
my_polygons_cells
##  class       : SpatVector 
##  geometry    : points 
##  dimensions  : 465, 4  (geometries, attributes)
##  extent      : 1, 50, -43, -1  (xmin, xmax, ymin, ymax)
##  coord. ref. :  
##  names       : sdimx sdimy cell_ID   poly_ID
##  type        : <int> <int>   <chr>     <chr>
##  values      :    10   -14   10x14 polygon 2
##                   10   -15   10x15 polygon 2
##                   10   -16   10x16 polygon 2

5.6 Differential Expression between selected regions

Find marker features for each polygon

scran_results <- findMarkers_one_vs_all(my_giotto_object,
                                        spat_unit = "cell",
                                        feat_type = "rna",
                                        method = "scran",
                                        expression_values = "normalized",
                                        cluster_column = "selections",
                                        min_feats = 10)
## using 'Scran' to detect marker feats. If used in published research, please cite:
##   Lun ATL, McCarthy DJ, Marioni JC (2016).
##   'A step-by-step workflow for low-level analysis of single-cell RNA-seq data with Bioconductor.'
##   F1000Res., 5, 2122. doi: 10.12688/f1000research.9501.2.
## 
##  start with cluster  no_polygon 
## 
##  start with cluster  polygon 1 
## 
##  start with cluster  polygon 2
top_genes <- scran_results[, head(.SD, 2), by = 'cluster']$feats

comparePolygonExpression will create a heatmap to observe the differential expression between polygon areas

comparePolygonExpression(my_giotto_object,
                         selected_feats = "top_genes")
## using 'Scran' to detect marker feats. If used in published research, please cite:
##   Lun ATL, McCarthy DJ, Marioni JC (2016).
##   'A step-by-step workflow for low-level analysis of single-cell RNA-seq data with Bioconductor.'
##   F1000Res., 5, 2122. doi: 10.12688/f1000research.9501.2.
## 
##  start with cluster  no_polygon 
## 
##  start with cluster  polygon 1 
## 
##  start with cluster  polygon 2

You can isolate each region to plot cell metadata such as the number of features per cell.

spatPlot2D(my_giotto_object,
           cell_color = 'nr_feats',
           group_by = 'selections',
           color_as_factor = FALSE,
           point_size = 2,
           coord_fix_ratio = 1,
           cow_n_col = 2,
           show_legend = TRUE)

Or plot gene expression per cell by selected region (polygon).

spatFeatPlot2D(my_giotto_object, 
               expression_values = 'scaled',
               feats = c("APOC1", "TMEM132D"), 
               group_by = 'selections',
               point_size = 2)

5.7 Cell types enriched between selected regions

You can compare the percent of cell types by selected region (polygon).

compareCellAbundance(my_giotto_object,
                     spat_unit = "cell",
                     feat_type = "rna")

You can isolate each region and plot cell types

##       cell_ID selections
##    1:   10x14  polygon 2
##    2:   10x15  polygon 2
##    3:   10x16  polygon 2
##    4:   10x17  polygon 2
##    5:   10x18 no_polygon
##   ---                   
## 1687:    9x46 no_polygon
## 1688:    9x47 no_polygon
## 1689:    9x48 no_polygon
## 1690:    9x49 no_polygon
## 1691:    9x50 no_polygon
spatPlot2D(my_giotto_object,
           spat_unit = "cell",
           feat_type = "protein",
           cell_color = 'integrated_leiden_clus',
           group_by = 'selections',
           color_as_factor = TRUE,
           point_size = 1,
           coord_fix_ratio = 1,
           cow_n_col = 2,
           show_legend = TRUE)

5.8 Plot polygons again

Finally, you can plot the regions previously selected.

plotPolygons(my_giotto_object,
             x = my_spatPlot,
             polygon_name = 'selections')